  MODULE TSTORMS_MOD
  implicit none

!=====================================================================
! --- NAMELIST
!=====================================================================

! integer, parameter :: gnx0 = 3   !for 1deg
 integer, parameter :: gnx0 = 6   !for 0.5deg
! integer, parameter :: gnx0 = 12  !for 0.25deg

 !gnx = 20 (winter storms)

  real :: crit_vort   =   5.0E-5   
  real :: crit_twc    =   0.8   
  real :: crit_thick  =  50.0   
  real :: crit_dist   =   4.0   
  real :: dist_twc    =   5.0
  real :: dist_thick  =   8.0
  real :: crit_psl    = 400.0
  real :: dist_psl    =   5.0
  real :: lat_bound_n =  90.0
  real :: lat_bound_s = -90.0
  real :: resolution  =   0.5
  real :: pslVortMaxDist = 8.0
  
  ! Added for GFDL tracker
  logical :: do_spline = .false.
  logical :: do_thickness = .false.
  
  ! pslVortMaxDist --- distance in degrees PSL and VORT are seperated by
  
  ! dist_psl --- distance in degrees for PSL falloff
  ! crit_psl --- amount of PSL (in Pa) rise from min
  ! crit_vort --- critical vorticity
  
  ! crit_dist --- how far away TWC and THICK need to be from PSL (degrees)
  ! dist_twc --- how far (in degrees) fall off for TWC needs to occur
  ! dist_thick --- how far (in degrees) fall off for THICK needs to occur
  
  ! crit_twc --- amount of TWC fall off in dist_twc is necessary (in K)
  ! crit_thick --- amount of THICK fall off in dist_thick is neces (in m)
  
  ! resolution --- grid resolution of model data (on uniform lat-lon grid)
  
  ! lat_bound_n, lat_bound_s --- where to not track vortices beyond

  namelist / nml_tstorms / crit_vort, crit_twc, crit_thick, crit_dist, &
                           dist_twc, crit_psl, dist_psl, dist_thick,   &
                           lat_bound_n, lat_bound_s, resolution, &
                           pslVortMaxDist

!=====================================================================
  contains

!######################################################################

  SUBROUTINE SET_TSTORMS
  implicit none
  READ( *, nml_tstorms )
  WRITE( *, nml_tstorms )
  end SUBROUTINE SET_TSTORMS

!######################################################################

!######################################################################

  SUBROUTINE TSTORMS ( Gwind, Gvort, Gtbar, Gpsl, Gthick,      &
                       Grlon, Grlat, iyear, imon, iday, ihour, &
                       iucy   )

!===================================================================
! --- LOCATE TROPICAL STORMS 
!===================================================================

  use TSGPAD_MOD, only : GPAD2,  GPAD1
  use SPLINE_MOD, only : SPLIE2, SPLIE3, FRPRM, SHAPE
  implicit none

!-------------------------------------------------------------------
! --- INPUT ARGUMENTS 
!     Gwind  - wind speed at 850 mb
!     Gvort  - vorticity  at 850 mb
!     Gtbar  - mean temperature for warm core layer
!     Gpsl   - sea level pressure
!     Gthick - thickness of 200 to 1000 mb layer 
!     Grlon  - longitudes
!     Grlat  - latitudes
!     iyear  - year
!     imon   - month
!     iday   - day of month
!     ihour  - hour
!     iucy   - unit for output
!-------------------------------------------------------------------
! --- OUTPUT - file "cyclones" 
!-------------------------------------------------------------------
! --- record # 1 
!     num0   - day
!     imon0  - month
!     iyear  - year
!     number - number of cyclones found
! --- records # 2...number+1 
!     idex, jdex - (i,j) index of cyclone 
!     svort_max  - max vorticity                  
!     swind_max  - max wind              
!      spsl_min  - min sea level pressure                 
!     svort_lon,  svort_lat - longitude & latitude of max vorticity 
!      spsl_lon,   spsl_lat - longitude & latitude of min slp 
!      stwc_lon,   stwc_lat - longitude & latitude of warm core 
!    sthick_lon, sthick_lat - longitude & latitude of max thickness 
!-------------------------------------------------------------------

  real,    intent(in),    dimension(:,:) :: Gwind, Gvort, Gtbar
  real,    intent(in),    dimension(:,:) :: Gpsl,  Gthick
  real,    intent(in),    dimension(:)   :: Grlon, Grlat
  integer, intent(in)                    :: ihour, iday,  imon, iyear, iucy

!-------------------------------------------------------------------
! --- LOCAL
!-------------------------------------------------------------------

  integer, parameter :: nx   = gnx0
! integer, parameter :: nx   = 12
  integer, parameter :: nx2  = 2*nx
  integer, parameter :: nxp1 = nx + 1

  real,    parameter :: ftol  = 0.01   
  integer, parameter :: nsmax = 10000

  real, dimension(SIZE(Grlon)+nx2) :: rlon
  real, dimension(SIZE(Grlat)+nx2) :: rlat

  real, dimension(SIZE(Gwind,1)+nx2,SIZE(Gwind,2)+nx2) ::  &
        vort, wind, psl,    tbar,    thick,                &
                    psl_dx, tbar_dx, thick_dx,             &
                    psl_dy, tbar_dy, thick_dy 

  real    :: vort_max, wind_max, psl_min, twc_max,   thick_max
  real    :: lon_vort,           lon_psl, lon_twc,   lon_thick
  real    :: lat_vort,           lat_psl, lat_twc,   lat_thick            
  logical ::                              exist_twc, exist_thick

  integer, dimension(nsmax) :: idex, jdex

  real,    dimension(nsmax) ::                                    &
           svort_max, swind_max,  spsl_min, stwc_max, sthick_max, &
           svort_lon,             spsl_lon, stwc_lon, sthick_lon, &
           svort_lat,             spsl_lat, stwc_lat, sthick_lat
  logical, dimension(nsmax) ::                                    &
                                            stwc_is,  sthick_is

  real,    dimension(2)  :: p

  real :: xx, yy, rr, fret

  integer :: ierr_pos, ierr_mag

  integer :: i, im, ip, ix, ixp3, ixp6
  integer :: j, jm, jp, jx, jxp3, jxp6, jxp6h 
  integer :: number, iter
  integer :: imm, jmm
  real    :: avg

!===================================================================

  ix    = SIZE( Gwind, 1 )
  jx    = SIZE( Gwind, 2 )
  ixp3  = ix + nx
  jxp3  = jx + nx
  ixp6  = ix + nx2
  jxp6  = jx + nx2
  jxp6h = jxp6 / 2

  number = 0

!-------------------------------------------------------------------
! --- SETUP
!-------------------------------------------------------------------

  CALL GPAD2( Gwind,  wind    )    !  Wind speed at 850 mb 
  CALL GPAD2( Gvort,  vort    )    !  Vorticity  at 850 mb 
  CALL GPAD2( Gtbar,  tbar    )    !  Mean temp for warm core layer
  CALL GPAD2( Gpsl,   psl     )    !  Sea level pressure.
if (do_thickness) then
  CALL GPAD2( Gthick, thick   )    !  Thickness of 200 to 1000 mb layer 
end if
!CMZ  CALL GPAD1( Grlon,  rlon, 0 )    !  Longitudes
!CMZ  CALL GPAD1( Grlat,  rlat, 1 )    !  Latitudes
  CALL GPAD1( Grlon,  rlon)    !  Longitudes
  CALL GPAD1( Grlat,  rlat)    !  Latitudes
 
! --- change sign of vorticity in southern hemisphere
  vort(:,1:jxp6h) = -1.0 * vort(:,1:jxp6h) 

! --- change sign of temperature & thickness
if (do_spline) then
   tbar(:,:) = -1.0 *  tbar(:,:)
if (do_thickness) then
  thick(:,:) = -1.0 * thick(:,:)
end if
end if

!-------------------------------------------------------------------
! --- INITIALIZE SPLINES
!-------------------------------------------------------------------
if (do_spline) then
  CALL SPLIE2( rlon, rlat, psl,   psl_dy   )
  CALL SPLIE3( rlon, rlat, psl,   psl_dx   )
 
  CALL SPLIE2( rlon, rlat, tbar,  tbar_dy  )
  CALL SPLIE3( rlon, rlat, tbar,  tbar_dx  )
 
  CALL SPLIE2( rlon, rlat, thick, thick_dy )
  CALL SPLIE3( rlon, rlat, thick, thick_dx )
end if

!===================================================================
! -- LOOP OVER GRID & LOOK FOR STORMS
!===================================================================

  do j = nxp1,jxp3
         if( ( rlat(j) > lat_bound_n ) .or. &
             ( rlat(j) < lat_bound_s ) ) CYCLE
  do i = nxp1,ixp3
! zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz

    im = i - nx
    ip = i + nx
    jm = j - nx
    jp = j + nx

!-------------------------------------------------------------------
! --- STEP 1: CHECK FOR VORTICITY MAX
!-------------------------------------------------------------------

   vort_max = MAXVAL( vort(im:ip,jm:jp) )
   wind_max = MAXVAL( wind(im:ip,jm:jp) )

! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  if( ( vort(i,j) /= vort_max  )  .or. &
      ( vort(i,j) <  crit_vort ) ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

  lon_vort = rlon(i)
  lat_vort = rlat(j)

!-------------------------------------------------------------------
! --- STEP 2: LOCATE LOCAL SEA LEVEL PRESSURE MIN 
!-------------------------------------------------------------------

  ierr_pos  = 0

  p(1) = lon_vort
  p(2) = lat_vort

if (do_spline) then
  CALL FRPRM( rlon, rlat, psl,  psl_dy, psl_dx, &
              p,    ftol, iter, fret,   ierr_pos )

  psl_min = fret
  lon_psl = p(1)
  lat_psl = p(2)
else
  call find_minmax (.true., psl(im:ip,jm:jp), imm, jmm, avg)
  imm=i-(nx+1)+imm
  jmm=j-(nx+1)+jmm
  lon_psl = rlon(imm)
  lat_psl = rlat(jmm)
  psl_min = psl (imm, jmm)
end if
  xx      = lon_psl - lon_vort
  yy      = lat_psl - lat_vort
  rr      = xx*xx + yy*yy

  if( rr >= crit_dist ) ierr_pos = 1

! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
   if ( ierr_pos == 1 ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

!-------------------------------------------------------------------
! --- STEP 3: CHECK FOR PRESENCE OF A WARM CORE
!-------------------------------------------------------------------

! --- location

  ierr_pos = 0
if (do_spline) then
  CALL FRPRM( rlon, rlat, tbar, tbar_dy, tbar_dx, &
              p,    ftol, iter, fret,    ierr_pos )

  if( ierr_pos == 0 ) then
    twc_max = -fret
    lon_twc = p(1)
    lat_twc = p(2)
    xx      = lon_twc - lon_psl
    yy      = lat_twc - lat_psl
    rr      = xx*xx + yy*yy
    if( rr >= crit_dist ) ierr_pos = 1
  else
    twc_max = 0.0
    lon_twc = 0.0
    lat_twc = 0.0
  endif

  exist_twc = (ierr_pos == 0)

! --- magnitude

  if ( exist_twc ) then

  ierr_mag = 0

  !CALL SHAPE( rlon, rlat, tbar, tbar_dy, tbar_dx,           &
  !            p,    ftol, iter, fret,    ierr_mag, crit_twc )
 
  exist_twc = exist_twc .and. ( ierr_mag == 0 )
 
  endif
else

  call find_minmax (.false., tbar(im:ip,jm:jp), imm, jmm, avg)
  imm=i-(nx+1)+imm
  jmm=j-(nx+1)+jmm
  lon_twc = rlon(imm)
  lat_twc = rlat(jmm)
  twc_max = tbar(imm, jmm) - avg
  xx      = lon_twc - lon_psl
  yy      = lat_twc - lat_psl
  rr      = xx*xx + yy*yy
  if( rr >= crit_dist ) ierr_pos = 1

  exist_twc = (ierr_pos == 0)

  exist_twc = exist_twc!miz .and. (twc_max > crit_twc)
!new
  exist_twc = exist_twc .and. (twc_max > crit_twc)
  if ( .not.exist_twc ) CYCLE
!new
end if

!-------------------------------------------------------------------
! --- STEP 4: CHECK FOR THICKNESS MAX
!-------------------------------------------------------------------

! --- location
if (do_thickness) then
  ierr_pos = 0

  CALL FRPRM( rlon, rlat, thick, thick_dy, thick_dx,  &
              p,    ftol, iter,  fret,     ierr_pos )

  if( ierr_pos == 0 ) then
    thick_max = -fret
    lon_thick = p(1)
    lat_thick = p(2)
    xx        = lon_thick - lon_psl
    yy        = lat_thick - lat_psl
    rr        = xx*xx + yy*yy
    if( rr >= crit_dist ) ierr_pos = 1
  else
    thick_max = 0.0
    lon_thick = 0.0
    lat_thick = 0.0
  endif

  exist_thick = (ierr_pos == 0)

! --- magnitude

  if ( exist_thick ) then

  ierr_mag = 0

  !CALL SHAPE( rlon, rlat, thick, thick_dy, thick_dx,     &
  !            p,    ftol, iter,  fret,     ierr_mag,  crit_thick )

  exist_thick = exist_thick .and. ( ierr_mag == 0 )
 
  endif
end if
!-------------------------------------------------------------------
! --- WE HAVE A TROPICAL STORM. SAVE INFO ABOUT STORM
!-------------------------------------------------------------------

    number = number + 1

  if( number > nsmax ) then
    PRINT *, '***************************************'
    PRINT *, '  GOT TOO MANY STORMS - INCREASE nsmax '
    PRINT *, '***************************************'
    STOP
  endif

        idex(number) = i - nx
        jdex(number) = j - nx

   svort_max(number) =  vort_max
   svort_lon(number) =  lon_vort
   svort_lat(number) =  lat_vort

   swind_max(number) =  wind_max

    spsl_min(number) =  psl_min
    spsl_lon(number) =  lon_psl
    spsl_lat(number) =  lat_psl

    stwc_max(number) =  twc_max
    stwc_lon(number) =  lon_twc
    stwc_lat(number) =  lat_twc
    stwc_is (number) =  exist_twc

if (.not.do_thickness) then
  sthick_max(number) =  10000
  sthick_lon(number) =  0.
  sthick_lat(number) =  0.
  sthick_is (number) =  .true.
end if

!miz  WRITE(*,*) i-nx, j-nx
! zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
  end do
  end do

!===================================================================
! ---  OUTPUT
!===================================================================

  WRITE(iucy,*) iday,  imon, iyear, number, ihour
  WRITE(   *,*) iyear, imon, iday,  ihour,  number

  if( number == 0 ) RETURN

  do i = 1,number
    WRITE(iucy,100) idex(i), jdex(i), spsl_lon(i), spsl_lat(i), &
              swind_max(i), svort_max(i), spsl_min(i), stwc_is(i), sthick_is(i)

!CMZ                stwc_is(i), sthick_is(i), stwc_max(i), sthick_max(i)
                
!10   FORMAT(I4, I4, F5.1, F5.1, F5.1, F5.1, F5.1)
100  FORMAT(I6,1x,I6,1x,2(f7.2,2x),F7.3,2x,ES10.3,2x,F10.2,2x,L1,1x,L1)

  end do

!===================================================================
  end SUBROUTINE TSTORMS

!######################################################################

!######################################################################

  SUBROUTINE TSTORMS_NOSPLINE ( Gwind, Gvort, Gtbar, Gpsl, Gthick, &
                       Grlon, Grlat, iyear, imon, iday, ihour, &
                       iucy )
                       
! This subroutine was added to deal with high resolution 0.25 deg
! or finer data that the spline routine takes too long to deal with
! Colin Zarzycki

!===================================================================
! --- LOCATE TROPICAL STORMS 
!===================================================================

  use TSGPAD_MOD, only : GPAD2,  GPAD1, GPAD1FLIP
  use SPLINE_MOD, only : SPLIE2, SPLIE3, FRPRM, SHAPE
  implicit none

!-------------------------------------------------------------------
! --- INPUT ARGUMENTS 
!     Gwind  - wind speed at 850 mb
!     Gvort  - vorticity  at 850 mb
!     Gtbar  - mean temperature for warm core layer
!     Gpsl   - sea level pressure
!     Gthick - thickness of 200 to 1000 mb layer 
!     Grlon  - longitudes
!     Grlat  - latitudes
!     iyear  - year
!     imon   - month
!     iday   - day of month
!     ihour  - hour
!     iucy   - unit for output
!-------------------------------------------------------------------
! --- OUTPUT - file "cyclones" 
!-------------------------------------------------------------------
! --- record # 1 
!     num0   - day
!     imon0  - month
!     iyear  - year
!     number - number of cyclones found
! --- records # 2...number+1 
!     idex, jdex - (i,j) index of cyclone 
!     svort_max  - max vorticity                  
!     swind_max  - max wind              
!      spsl_min  - min sea level pressure                 
!     svort_lon,  svort_lat - longitude & latitude of max vorticity 
!      spsl_lon,   spsl_lat - longitude & latitude of min slp 
!      stwc_lon,   stwc_lat - longitude & latitude of warm core 
!    sthick_lon, sthick_lat - longitude & latitude of max thickness 
!-------------------------------------------------------------------

  real,    intent(in),    dimension(:,:) :: Gwind, Gvort,  Gtbar
  real,    intent(in),    dimension(:,:) :: Gpsl,  Gthick
  real,    intent(in),    dimension(:)   :: Grlon, Grlat
  integer, intent(in)                    :: ihour, iday,  imon, iyear, iucy

!-------------------------------------------------------------------
! --- LOCAL
!-------------------------------------------------------------------

  integer, parameter :: nx   = gnx0
  !integer, parameter :: nx   = 16
  integer, parameter :: nx2  = 2*nx
  integer, parameter :: nxp1 = nx + 1

  real,    parameter :: ftol  = 0.01   
  integer, parameter :: nsmax = 10000

  real, dimension(SIZE(Grlon)+nx2) :: rlon
  real, dimension(SIZE(Grlat)+nx2) :: rlat

  real, dimension(SIZE(Gwind,1)+nx2,SIZE(Gwind,2)+nx2) ::  &
        vort, wind, psl,    tbar,    thick,            &
                    psl_dx, tbar_dx, thick_dx,         &
                    psl_dy, tbar_dy, thick_dy 

  real    :: vort_max, wind_max, psl_min, twc_max,   thick_max
  real    :: lon_vort,           lon_psl, lon_twc,   lon_thick
  real    :: lat_vort,           lat_psl, lat_twc,   lat_thick            
  logical ::                              exist_twc, exist_thick

  integer, dimension(nsmax) :: idex, jdex

  real,    dimension(nsmax) ::                                    &
           svort_max, swind_max,  spsl_min, stwc_max, sthick_max, &
           svort_lon,             spsl_lon, stwc_lon, sthick_lon, &
           svort_lat,             spsl_lat, stwc_lat, sthick_lat, &
           spsl_fail
           
  logical, dimension(nsmax) ::                                    &
                                            stwc_is,  sthick_is

  real,    dimension(2)  :: p

  real :: xx, yy, rr, fret

  integer :: ierr_pos, ierr_mag

  integer :: i, im, ip, ix, ixp3, ixp6
  integer :: j, jm, jp, jx, jxp3, jxp6, jxp6h 
  integer :: number, iter
  
!====== Colin stuff

  integer :: starti,startj
  
  real :: distbwpts,avg

  integer :: ii,jj
  
  integer :: psmin_i,psmin_j,psl_crit_dist,num_fail_psl
  integer :: psl_off_dist, psl_od_diag
  real    :: psmn

  integer :: tbrmx_i,tbrmx_j,twc_crit_dist,num_fail_twc
  integer :: twc_off_dist, twc_od_diag
  real    :: tbar_max
  
  integer :: tckmx_i,tckmx_j,thick_crit_dist,num_fail_thick
  integer :: thick_off_dist, thick_od_diag
  real    :: tck_max

!  integer, allocatable :: twcPSLColMtx(:,:), pslVortColMtx(:,:)

!===================================================================

  ix    = SIZE( Gwind, 1 )
  jx    = SIZE( Gwind, 2 )
  ixp3  = ix + nx
  jxp3  = jx + nx
  ixp6  = ix + nx2
  jxp6  = jx + nx2
  jxp6h = jxp6 / 2

  number = 0

!-------------------------------------------------------------------
! --- SETUP
!-------------------------------------------------------------------

! Here we are "padding" the arrays with ghost cells so we can do the
! vorticity and wind checks at all locations

  CALL GPAD2( Gwind,  wind  )    !  Wind speed at 850 mb 
  CALL GPAD2( Gvort,  vort  )    !  Vorticity  at 850 mb 
  CALL GPAD2( Gtbar,  tbar  )    !  Mean temp for warm core layer
  CALL GPAD2( Gpsl,   psl   )    !  Sea level pressure.
  CALL GPAD2( Gthick, thick )    !  Thickness of 200 to 1000 mb layer 
  !CALL GPAD1FLIP( Grlon,  rlon  , 0)    !  Longitudes
  !CALL GPAD1FLIP( Grlat,  rlat  , 1)    !  Latitudes  
  CALL GPAD1( Grlon,  rlon)    !  Longitudes
  CALL GPAD1( Grlat,  rlat)    !  Latitudes 
 
! --- change sign of vorticity in southern hemisphere
  do j=1,jxp6
  if( rlat(j) < 0.0 ) then
    vort(:,j) = -1.0 * vort(:,j) 
  end if
  end do

! --- change sign of temperature & thickness
!   tbar(:,:) = -1.0 *  tbar(:,:)
!  thick(:,:) = -1.0 * thick(:,:)

!===================================================================
! -- Define areas
!===================================================================
  
  print *, '---------------------'
  print *, 'resolution ', resolution
  print *, 'pslVortMaxDist ',pslVortMaxDist
  print *, 'crit_dist ', crit_dist
  print *, 'crit_vort ', crit_vort
  print *, 'crit_twc ', crit_twc
  print *, 'dist_twc ', dist_twc
  print *, 'crit_thick ', crit_thick
  print *, 'dist_thick ', dist_thick
  print *, 'crit_psl ', crit_psl
  print *, 'dist_psl ', dist_psl
  print *, '---------------------'
  
  psl_crit_dist = CEILING(pslVortMaxDist/resolution)
  print *, 'psl_crit_dist ', psl_crit_dist

  psl_off_dist = CEILING(dist_psl/resolution)
  psl_od_diag = CEILING(psl_off_dist/SQRT(2.))
  print *, 'psl_off_dist ', psl_off_dist
  print *, 'psl_od_diag ', psl_od_diag

  twc_crit_dist = CEILING(crit_dist/resolution)
  print *, 'twc_crit_dist', twc_crit_dist
  twc_off_dist = CEILING(dist_twc/resolution)
  twc_od_diag = CEILING(twc_off_dist/SQRT(2.))
  print *, 'twc_off_dist ', twc_off_dist
  print *, 'twc_od_diag ', twc_od_diag
  
  thick_crit_dist = CEILING(crit_dist/resolution)
  print *, 'thick_crit_dist', thick_crit_dist
  thick_off_dist = CEILING(dist_thick/resolution)
  thick_od_diag = CEILING(thick_off_dist/SQRT(2.))
  print *, 'thick_off_dist ', thick_off_dist
  print *, 'thick_od_diag ', thick_od_diag
  
!-------------------------------------------------------------------
! --- CREATE TWC COLLOCATION MATRICES (COLIN)
!-------------------------------------------------------------------
  
!  allocate ( twcPSLColMtx(2*twc_crit_dist+1,2*twc_crit_dist+1) )
!  call circleMatrix(twcPSLColMtx, twc_crit_dist)

!-------------------------------------------------------------------
! --- CREATE PSL-VORT DISTANCE MATRICES
!-------------------------------------------------------------------

! This one is currently used for PSL - vort distance
!  allocate ( pslVortColMtx(2*psl_crit_dist+1,2*psl_crit_dist+1) )
!  call circleMatrix(pslVortColMtx, psl_crit_dist)
  
!===================================================================
! -- LOOP OVER GRID & LOOK FOR STORMS
!===================================================================

  print *, 'nx ', nx
  print *, 'nx2 ', nx2
  print *, 'nxp1 ', nxp1
  print *, 'ix ', ix
  print *, 'jx ', jx
  print *, 'ixp3 ', ixp3
  print *, 'jxp3 ', jxp3
  print *, 'ixp6 ', ixp6
  print *, 'jxp6 ', jxp6
  print *, 'jxp6h ', jxp6h
  
  do j = nxp1,jxp3
         if( ( rlat(j) > lat_bound_n ) .or. &
             ( rlat(j) < lat_bound_s ) ) CYCLE
  do i = nxp1,ixp3
! zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz

    im = i - nx
    ip = i + nx
    jm = j - nx
    jp = j + nx

!-------------------------------------------------------------------
! --- STEP 1: CHECK FOR VORTICITY MAX
!-------------------------------------------------------------------

   vort_max = MAXVAL( vort(im:ip,jm:jp) )
   wind_max = MAXVAL( wind(im:ip,jm:jp) )

! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  if( ( vort(i,j) /= vort_max  )  .or. &
      ( vort(i,j) <  crit_vort ) ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

  lon_vort = rlon(i)
  lat_vort = rlat(j)
  
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
! Make sure vorticity max is not at equator
!  if( ABS(lat_vort) < 1.1  ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  
  PRINT *, '        '
  PRINT *, '----------------------'
  PRINT *, 'O.O.O.O - FOUND VORT MAX at ',lon_vort,' ',lat_vort

!-------------------------------------------------------------------
! --- STEP 2: LOCATE LOCAL SEA LEVEL PRESSURE MIN 
!-------------------------------------------------------------------

  ierr_pos  = 0
  
! Set large number for initial psmn
  psmn = 1E7
  psmin_i = 1
  psmin_j = 1
 
! Find starting indices for loop (top left corner of square)
  starti = i-psl_crit_dist
  startj = j-psl_crit_dist
  
!- ! Loop over gridboxes within critdist to see where psl_min is
!-   do ii = i-psl_crit_dist,i+psl_crit_dist
!-   !write (*,'(<33>(I1,X))') ( pslVortColMtx(ii-starti+1,jj-startj+1), jj=j-psl_crit_dist,j+psl_crit_dist )
!-   do jj = j-psl_crit_dist,j+psl_crit_dist
!-   !Makes sure it's within circular domain defined by pslVortColMtx
!-   !print *,ii-starti+1,' ',jj-startj+1
!-  !if (pslVortColMtx(ii-starti+1,jj-startj+1) .gt. 0) then
!-   if( SQRT((ii-i)**2.+(jj-j)**2.) .le. psl_crit_dist) then
!-       if (psl(ii,jj).lt.psmn) then
!-         psmn = psl(ii,jj)
!-         psmin_i = ii
!-         psmin_j = jj
!-       else
!-       end if
!-   end if
!-   end do
!-   end do

  call find_minmax_gen (.true.,psl(i-psl_crit_dist:i+psl_crit_dist,j-psl_crit_dist:j+psl_crit_dist),ii,jj,avg)
  psmin_i = ii+starti-1
  psmin_j = jj+startj-1
  psmn = psl(psmin_i,psmin_j)

  lon_psl = rlon(psmin_i)
  lat_psl = rlat(psmin_j)
  print *, '  FOUND PSL MIN at ',lon_psl,' ',lat_psl

  ! Calculate distance between vort center (i,j) and psl (psmin_i,psmin_j)
  ! in grid boxes
  distbwpts = SQRT( (psmin_i-i)**2. + (psmin_j-j)**2. )
  if (distbwpts .gt. psl_crit_dist) then
    print *, 'X.X.X.X - PSL dist = ',distbwpts,' compared to ',psl_crit_dist
    ierr_pos = 1
  else
    print *, 'O.O.O.O - PSL dist = ',distbwpts,' compared to ',psl_crit_dist
    ierr_pos = 0
  end if
  
  psl_min = psmn
  
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
   !if ( ierr_pos == 1 ) print *, 'TOSSING PSL POS', lon_psl, lat_psl
   if ( ierr_pos == 1 ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

  num_fail_psl=0
  ierr_mag=0
  call check_contour(.true.,psl,psmin_i,psmin_j,num_fail_psl,0,ierr_mag,psl_off_dist,psl_od_diag,crit_psl)
  
  if ( ierr_mag == 1 ) then
    print *, 'X.X.X.X - PSL NUMFAIL = ',num_fail_psl,' compared to 0'
  else
    print *, 'O.O.O.O - PSL NUMFAIL = ',num_fail_psl,' compared to 0'
  end if

! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
   !if ( ierr_mag == 1 ) print *, 'TOSSING  PSL MAG', lon_psl, lat_psl, ' BECAUSE NUMFAIL ',num_fail_psl
   if ( ierr_mag == 1 ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

!-------------------------------------------------------------------
! --- STEP 3: CHECK FOR PRESENCE OF A WARM CORE
!-------------------------------------------------------------------

  ierr_pos = 0
  
  tbar_max = 0
  tbrmx_i = 1
  tbrmx_j = 1
  
  starti = psmin_i-twc_crit_dist
  startj = psmin_j-twc_crit_dist
  
  do ii = psmin_i-twc_crit_dist,psmin_i+twc_crit_dist
  do jj = psmin_j-twc_crit_dist,psmin_j+twc_crit_dist
 !  if (twcPSLColMtx(ii-starti+1,jj-startj+1) .gt. 0) then
  if( SQRT((psmin_i-ii)**2.+(psmin_j-jj)**2.) .le. twc_crit_dist) then
      if (tbar(ii,jj).gt.tbar_max) then
        tbar_max = tbar(ii,jj)
        tbrmx_i = ii
        tbrmx_j = jj
      else
      end if
  end if
  end do
  end do

!  call find_minmax_gen (.false., tbar(psmin_i-twc_crit_dist:psmin_i+twc_crit_dist,psmin_j-twc_crit_dist:psmin_j+twc_crit_dist), ii, jj)
!  tbrmx_i = ii+starti-1
!  tbrmx_j = jj+startj-1
!  tbar_max = tbar(tbrmx_i,tbrmx_j)

  lon_twc = rlon(tbrmx_i)
  lat_twc = rlat(tbrmx_j)
  print *, '  FOUND TWC MAX at ',lon_twc,' ',lat_twc
  
  distbwpts = SQRT( (psmin_i-tbrmx_i)**2. + (psmin_j-tbrmx_j)**2. )

  if (distbwpts .gt. twc_crit_dist) then
    print *, 'X.X.X.X - TWC dist = ',distbwpts,' compared to ',twc_crit_dist
    ierr_pos = 1
  else
    print *, 'O.O.O.O - TWC dist = ',distbwpts,' compared to ',twc_crit_dist
    ierr_pos = 0
  end if
  
  twc_max = tbar_max

  exist_twc = ( ierr_pos == 0 )

! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
   !if ( ierr_pos == 1 ) print *, 'TOSSING TBAR POS', lon_psl, lat_psl
   if ( ierr_pos == 1 ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

! --- magnitude

  num_fail_twc=0
  ierr_mag=0
  call check_contour(.false.,tbar,tbrmx_i,tbrmx_j,num_fail_twc,0,ierr_mag,twc_off_dist,twc_od_diag,crit_twc)

  if ( ierr_mag == 1 ) then
    print *, 'X.X.X.X - TWC NUMFAIL = ',num_fail_twc,' compared to 1'
  else
    print *, 'O.O.O.O - TWC NUMFAIL = ',num_fail_twc,' compared to 1'
  end if

  exist_twc = exist_twc .and. ( ierr_mag == 0 )

! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
   !if ( ierr_mag == 1 ) print *, 'TOSSING  TBAR MAG', lon_psl, lat_psl, ' BECAUSE NUMFAIL ',num_fail_twc
   if ( ierr_mag == 1 ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

!-------------------------------------------------------------------
! --- WE HAVE A TROPICAL STORM. SAVE INFO ABOUT STORM
!-------------------------------------------------------------------

    number = number + 1
    
    print *, 'ooooooooooooo -- CANDIDATE CYCLONE'
    print *, number, lon_psl, lat_psl!, '     ', lon_twc, lat_twc
    print *, ' '

  if( number > nsmax ) then
    PRINT *, '***************************************'
    PRINT *, '  GOT TOO MANY STORMS - INCREASE nsmax '
    PRINT *, '***************************************'
    STOP
  endif

        idex(number) = i - nx
        jdex(number) = j - nx

   svort_max(number) =  vort_max
   svort_lon(number) =  lon_vort
   svort_lat(number) =  lat_vort

   swind_max(number) =  wind_max

    spsl_min(number) =  psl_min
    spsl_lon(number) =  lon_psl
    spsl_lat(number) =  lat_psl
   spsl_fail(number) =  num_fail_psl

    stwc_max(number) =  twc_max
    stwc_lon(number) =  lon_twc
    stwc_lat(number) =  lat_twc
    stwc_is (number) =  exist_twc

  sthick_max(number) =  thick_max
  sthick_lon(number) =  lon_thick
  sthick_lat(number) =  lat_thick
  sthick_is (number) =  exist_thick

!  WRITE(*,*) i-nx, j-nx
! zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
  end do
  end do

!===================================================================
! ---  OUTPUT
!===================================================================

  WRITE(iucy,*) iday,  imon, iyear, number, ihour
  WRITE(   *,*) iyear, imon, iday,  ihour,  number

  if( number == 0 ) RETURN

! ---  sthick_is(1:number) = .true.

!  do i = 1,number
!    WRITE(iucy,*) idex(i),       jdex(i),              & 
!              spsl_lon(i),   spsl_lat(i),              &
!              swind_max(i), svort_max(i), spsl_min(i), &
!                stwc_is(i), sthick_is(i)
!CMZ                stwc_is(i), sthick_is(i), spsl_fail(i), spsl_drop(i)

  do i = 1,number
    WRITE(iucy,100) idex(i), jdex(i), spsl_lon(i), spsl_lat(i), &
              swind_max(i), svort_max(i), spsl_min(i), stwc_is(i), sthick_is(i)

!CMZ                stwc_is(i), sthick_is(i), stwc_max(i), sthick_max(i)
                
!10   FORMAT(I4, I4, F5.1, F5.1, F5.1, F5.1, F5.1)
100  FORMAT(I6,1x,I6,1x,2(f7.2,2x),F7.3,2x,ES10.3,2x,F10.2,2x,L1,1x,L1)

  end do
  
!  deallocate(twcPSLColMtx)
!  deallocate(pslVortColMtx)
  
!===================================================================
  end SUBROUTINE TSTORMS_NOSPLINE
  
!######################################################################

!######################################################################
  
  SUBROUTINE TSTORMS_HODGES ( Gwind, Gvort, Gvort850T42, Gvort200, &
                   Gvort500, &
                   Gpsl, Grlon, Grlat, iyear, imon, iday, ihour, &
                    iucy )

! This subroutine was added to do Hodges tracker

!===================================================================
! --- LOCATE TROPICAL STORMS 
!===================================================================

  use TSGPAD_MOD, only : GPAD2,  GPAD1
  use SPLINE_MOD, only : SPLIE2, SPLIE3, FRPRM, SHAPE
  implicit none

!-------------------------------------------------------------------
! --- INPUT ARGUMENTS 
!     Gwind  - wind speed at 850 mb
!     Gvort  - vorticity  at 850 mb
!     Gtbar  - mean temperature for warm core layer
!     Gpsl   - sea level pressure
!     Gthick - thickness of 200 to 1000 mb layer 
!     Grlon  - longitudes
!     Grlat  - latitudes
!     iyear  - year
!     imon   - month
!     iday   - day of month
!     ihour  - hour
!     iucy   - unit for output
!-------------------------------------------------------------------
! --- OUTPUT - file "cyclones" 
!-------------------------------------------------------------------
! --- record # 1 
!     num0   - day
!     imon0  - month
!     iyear  - year
!     number - number of cyclones found
! --- records # 2...number+1 
!     idex, jdex - (i,j) index of cyclone 
!     svort_max  - max vorticity                  
!     swind_max  - max wind              
!      spsl_min  - min sea level pressure                 
!     svort_lon,  svort_lat - longitude & latitude of max vorticity 
!      spsl_lon,   spsl_lat - longitude & latitude of min slp 
!      stwc_lon,   stwc_lat - longitude & latitude of warm core 
!    sthick_lon, sthick_lat - longitude & latitude of max thickness 
!-------------------------------------------------------------------

  real,    intent(in),    dimension(:,:) :: Gwind, Gvort, Gvort200, Gvort500
  real,    intent(in),    dimension(:,:) :: Gvort850T42
  real,    intent(in),    dimension(:,:) :: Gpsl
  real,    intent(in),    dimension(:)   :: Grlon, Grlat
  integer, intent(in)                    :: ihour, iday,  imon, iyear, iucy

!-------------------------------------------------------------------
! --- LOCAL
!-------------------------------------------------------------------

  integer, parameter :: nx   = 16
  integer, parameter :: nx2  = 2*nx
  integer, parameter :: nxp1 = nx + 1

  real,    parameter :: ftol  = 0.01   
  integer, parameter :: nsmax = 10000

  real, dimension(SIZE(Grlon)+nx2) :: rlon
  real, dimension(SIZE(Grlat)+nx2) :: rlat

  real, dimension(SIZE(Gwind,1)+nx2,SIZE(Gwind,2)+nx2) ::  &
        vort, wind, psl, vort850T42, vort500,  vort200,     &
                    psl_dx, tbar_dx, thick_dx,             &
                    psl_dy, tbar_dy, thick_dy 

  real    :: vort_max, wind_max, psl_min, twc_max,   thick_max
  real    :: lon_vort,           lon_psl, lon_twc,   lon_thick
  real    :: lat_vort,           lat_psl, lat_twc,   lat_thick            
  logical ::                              exist_twc, exist_thick

  integer, dimension(nsmax) :: idex, jdex

  real,    dimension(nsmax) ::                                    &
           svort_max, swind_max,  spsl_min, stwc_max, sthick_max, &
           svort_lon,             spsl_lon, stwc_lon, sthick_lon, &
           svort_lat,             spsl_lat, stwc_lat, sthick_lat, &
           spsl_fail, spsl_drop
           
  logical, dimension(nsmax) ::                                    &
                                            stwc_is,  sthick_is

  real,    dimension(2)  :: p

  real :: xx, yy, rr, fret

  integer :: ierr_pos, ierr_mag

  integer :: i, im, ip, ix, ixp3, ixp6
  integer :: j, jm, jp, jx, jxp3, jxp6, jxp6h 
  integer :: number, iter
  
!====== Colin stuff

  integer :: starti,startj
  
  real :: distbwpts
  
  integer :: ii,jj,psmin_i,psmin_j,psl_crit_dist
  integer :: vr200mx_i,vr200mx_j
  integer :: vr500mx_i,vr500mx_j
  integer :: vr850mx_i,vr850mx_j
  real    :: psmn, vr200mx, vr500mx, vr850mx
  real    :: vort200_max, vort500_max, vort850_max
  real    :: vr_diff, min_vr_diff, vr_sign
  real    :: lon_vr200mx, lat_vr200mx
  real    :: lon_vr500mx, lat_vr500mx
  real    :: lon_vr850mx, lat_vr850mx
  real    :: vr_crit_dist

  
!  real :: resolution, eps, dist
  real :: eps, dist
  integer :: degree_thresh, thresh_grid_dist, center_i, center_j
  integer :: iii, jjj, count
  integer, allocatable :: twcPSLColMtx(:,:), pslVortColMtx(:,:), ringmtx(:,:)

!===================================================================

  ix    = SIZE( Gwind, 1 )
  jx    = SIZE( Gwind, 2 )
  ixp3  = ix + nx
  jxp3  = jx + nx
  ixp6  = ix + nx2
  jxp6  = jx + nx2
  jxp6h = jxp6 / 2

  number = 0

!-------------------------------------------------------------------
! --- SETUP
!-------------------------------------------------------------------

! Here we are "padding" the arrays with ghost cells so we can do the
! vorticity and wind checks at all locations

  CALL GPAD2( Gwind,  wind  )        !  Wind speed at 850 mb 
  CALL GPAD2( Gvort,  vort  )        !  Vorticity  at 850 mb
  CALL GPAD2( Gvort850T42,  vort850T42  )  !  Vorticity  at 200 mb 
  CALL GPAD2( Gvort200,  vort200  )  !  Vorticity  at 200 mb
  CALL GPAD2( Gvort500,  vort500  )  !  Vorticity  at 200 mb 
  CALL GPAD2( Gpsl,   psl   )        !  Sea level pressure.
  CALL GPAD1( Grlon,  rlon  )        !  Longitudes
  CALL GPAD1( Grlat,  rlat  )        !  Latitudes  
 
! --- change sign of vorticity in southern hemisphere
  do j=1,jxp6
  if( rlat(j) < 0.0 ) then
    vort(:,j)    = -1.0 * vort(:,j)
    vort850T42(:,j) = -1.0 * vort850T42(:,j) 
    vort200(:,j) = -1.0 * vort200(:,j) 
    vort500(:,j) = -1.0 * vort500(:,j) 
  end if
  end do

! --- change sign of temperature & thickness
!   tbar(:,:) = -1.0 *  tbar(:,:)
!  thick(:,:) = -1.0 * thick(:,:)

  eps = 10e-15

!===================================================================
! -- Define areas
!===================================================================
  
  print *, '---------------------'
  print *, 'resolution ', resolution
  print *, 'pslVortMaxDist ',pslVortMaxDist
  print *, 'crit_dist ', crit_dist
  print *, 'crit_vort ', crit_vort
  print *, 'crit_twc ', crit_twc
  print *, 'dist_twc ', dist_twc
  print *, 'crit_thick ', crit_thick
  print *, 'dist_thick ', dist_thick
  print *, 'crit_psl ', crit_psl
  print *, 'dist_psl ', dist_psl
  print *, '---------------------'
  

  
  psl_crit_dist = CEILING(pslVortMaxDist/resolution)
  print *, 'psl_crit_dist ', psl_crit_dist
  
  vr_crit_dist = psl_crit_dist
  
!-------------------------------------------------------------------
! --- CREATE PSL-VORT DISTANCE MATRICES
!-------------------------------------------------------------------

! This one is currently used for PSL - vort distance

  thresh_grid_dist = psl_crit_dist

  center_i = thresh_grid_dist+1
  center_j = thresh_grid_dist+1

  allocate ( pslVortColMtx(2*thresh_grid_dist+1,2*thresh_grid_dist+1) )
  
  do iii = 1,2*thresh_grid_dist+1
    do jjj = 1,2*thresh_grid_dist+1
      dist = ( (iii-center_i)**2 + (jjj-center_j)**2 )
      dist = dist ** (0.5)
      if (dist .LE. (thresh_grid_dist + eps)) then
        pslVortColMtx(iii,jjj) = 1
      else
        pslVortColMtx(iii,jjj) = 0
      end if
    end do
  end do
  
!  do iii=1,2*thresh_grid_dist+1
!    write (*,'(<2*thresh_grid_dist+1>(I1,X))') ( pslVortColMtx(iii,jjj), jjj=1,2*thresh_grid_dist+1 )
!  enddo
  
!===================================================================
! -- LOOP OVER GRID & LOOK FOR STORMS
!===================================================================

  print *, 'nx ', nx
  print *, 'nx2 ', nx2
  print *, 'nxp1 ', nxp1
  print *, 'ix ', ix
  print *, 'jx ', jx
  print *, 'ixp3 ', ixp3
  print *, 'jxp3 ', jxp3
  print *, 'ixp6 ', ixp6
  print *, 'jxp6 ', jxp6
  print *, 'jxp6h ', jxp6h
  
  do j = nxp1,jxp3
         if( ( rlat(j) > lat_bound_n ) .or. &
             ( rlat(j) < lat_bound_s ) ) CYCLE
  do i = nxp1,ixp3
! zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz

    im = i - nx
    ip = i + nx
    jm = j - nx
    jp = j + nx

!-------------------------------------------------------------------
! --- STEP 1: CHECK FOR VORTICITY MAX
!-------------------------------------------------------------------

   vort_max = MAXVAL( vort850T42(im:ip,jm:jp) )
   wind_max = MAXVAL( wind(im:ip,jm:jp) )

! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  if( ( vort850T42(i,j) /= vort_max  )  .or. &
      ( vort850T42(i,j) <  crit_vort ) ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

  lon_vort = rlon(i)
  lat_vort = rlat(j)
  
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
! Make sure vorticity max is not at equator
  if( ABS(lat_vort) < 1.1  ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  
  PRINT *, '        '
  PRINT *, '----------------------'
  PRINT *, 'O.O.O.O - FOUND 850 MB VORT MAX at ',lon_vort,' ',lat_vort

  
!-------------------------------------------------------------------
! --- STEP 2a: LOCATE 850 MB VORT MAX
!-------------------------------------------------------------------

  ierr_pos  = 0
  
! Set large number for initial psmn
  vr850mx = 0
  vr850mx_i = 1
  vr850mx_j = 1
 
! Find starting indices for loop (top left corner of square)
  starti = i-vr_crit_dist
  startj = j-vr_crit_dist
  
! Loop over gridboxes within critdist to see where psl_min is
  do ii = i-vr_crit_dist,i+vr_crit_dist
  do jj = j-vr_crit_dist,j+vr_crit_dist
      if (vort(ii,jj).gt.vr850mx) then
        vr850mx = vort(ii,jj)
        vr850mx_i = ii
        vr850mx_j = jj
      else
      end if
  end do
  end do
  
  lon_vr850mx = rlon(vr850mx_i)
  lat_vr850mx = rlat(vr850mx_j)
  
!  if (vr850mx .le. eps) then
!    print *, 'X.X.X.X - DID NOT FIND 850 MB VORT MAX'
!    CYCLE
!  else
!    PRINT *, 'O.O.O.O - FOUND 850 MB VORT MAX at ',lon_vr850mx,' ',lat_vr850mx
!  end if 
  
  ! Calculate distance between vort center (i,j) and psl (psmin_i,psmin_j)
  ! in grid boxes
  distbwpts = SQRT( (vr850mx_i-i)**2. + (vr850mx_j-j)**2. )
  
!  if (distbwpts .gt. vr_crit_dist) then
!    print *, 'X.X.X.X - PSL dist = ',distbwpts,' compared to ',psl_crit_dist
!    ierr_pos = 1
!  else
!    print *, 'O.O.O.O - PSL dist = ',distbwpts,' compared to ',psl_crit_dist
!    ierr_pos = 0
!  end if

  vort850_max = vr850mx  
  PRINT *, vort850_max  
  
  
  
  
  
!-------------------------------------------------------------------
! --- STEP 2a: LOCATE 500 MB VORT MAX
!-------------------------------------------------------------------

  ierr_pos  = 0
  
! Set large number for initial psmn
  vr500mx = 0
  vr500mx_i = 1
  vr500mx_j = 1
 
! Find starting indices for loop (top left corner of square)
  starti = i-vr_crit_dist
  startj = j-vr_crit_dist
  
! Loop over gridboxes within critdist to see where psl_min is
  do ii = i-vr_crit_dist,i+vr_crit_dist
  do jj = j-vr_crit_dist,j+vr_crit_dist
      if (vort500(ii,jj).gt.vr500mx) then
        vr500mx = vort500(ii,jj)
        vr500mx_i = ii
        vr500mx_j = jj
      else
      end if
  end do
  end do
  
  lon_vr500mx = rlon(vr500mx_i)
  lat_vr500mx = rlat(vr500mx_j)
  
  if (vr500mx .le. eps) then
    print *, 'X.X.X.X - DID NOT FIND 500 MB VORT MAX'
    CYCLE
  else
    PRINT *, 'O.O.O.O - FOUND 500 MB VORT MAX at ',lon_vr500mx,' ',lat_vr500mx
    !PRINT *, 'vr500mx ',vr500mx
  end if 
  
  ! Calculate distance between vort center (i,j) and psl (psmin_i,psmin_j)
  ! in grid boxes
  distbwpts = SQRT( (vr500mx_i-i)**2. + (vr500mx_j-j)**2. )
  
!  if (distbwpts .gt. vr_crit_dist) then
!    print *, 'X.X.X.X - PSL dist = ',distbwpts,' compared to ',psl_crit_dist
!    ierr_pos = 1
!  else
!    print *, 'O.O.O.O - PSL dist = ',distbwpts,' compared to ',psl_crit_dist
!    ierr_pos = 0
!  end if

  vort500_max = vr500mx
  
!-------------------------------------------------------------------
! --- STEP 2b: LOCATE 200 MB VORT MAX
!-------------------------------------------------------------------

  ierr_pos  = 0
  
! Set large number for initial psmn
  vr200mx = 0
  vr200mx_i = 1
  vr200mx_j = 1
 
! Find starting indices for loop (top left corner of square)
  starti = i-vr_crit_dist
  startj = j-vr_crit_dist
  
! Loop over gridboxes within critdist to see where psl_min is
  do ii = i-vr_crit_dist,i+vr_crit_dist
  do jj = j-vr_crit_dist,j+vr_crit_dist
      if (vort200(ii,jj).gt.vr200mx) then
        vr200mx = vort200(ii,jj)
        vr200mx_i = ii
        vr200mx_j = jj
      else
      end if
  end do
  end do
  
  lon_vr200mx = rlon(vr200mx_i)
  lat_vr200mx = rlat(vr200mx_j)
  
  if (vr200mx .le. eps) then
    print *, 'X.X.X.X - DID NOT FIND 200 MB VORT MAX'
    CYCLE
  else
    PRINT *, 'O.O.O.O - FOUND 200 MB VORT MAX at ',lon_vr200mx,' ',lat_vr200mx
    !PRINT *, 'vr200mx ',vr200mx
  end if 
  
  ! Calculate distance between vort center (i,j) and psl (psmin_i,psmin_j)
  ! in grid boxes
  distbwpts = SQRT( (vr200mx_i-i)**2. + (vr200mx_j-j)**2. )
  
!  if (distbwpts .gt. vr_crit_dist) then
!    print *, 'X.X.X.X - PSL dist = ',distbwpts,' compared to ',psl_crit_dist
!    ierr_pos = 1
!  else
!    print *, 'O.O.O.O - PSL dist = ',distbwpts,' compared to ',psl_crit_dist
!    ierr_pos = 0
!  end if

  vort200_max = vr200mx
  
!-------------------------------------------------------------------
! --- STEP 3: DOES MINIMUM REDUCTION OCCUR? 
!-------------------------------------------------------------------

  vr_diff = 0
  vr_diff = vort850_max - vort200_max
  min_vr_diff = 6.0e-5
  
  if (vr_diff .le. min_vr_diff) then
    print *, 'X.X.X.X - VORT DIFFERENCE TOO SMALL ',vr_diff,' < ',min_vr_diff
    CYCLE
  else
    PRINT *, 'O.O.O.O - FOUND VORT DIFF ',vr_diff,' > ',min_vr_diff
  end if   

!-------------------------------------------------------------------
! --- STEP 3: DOES VORTICITY DECREASE WITH HEIGHT?
!-------------------------------------------------------------------

  vr_sign = 0
  vr_sign = vort850_max - vort500_max
  
  if (vr_sign .le. 0) then
    print *, 'X.X.X.X - VORT INCREASES BW 850 and 500: ',vr_sign
    print *, 'X.X.X.X - 850 ',vort_max,'     500 ',vort500_max
    CYCLE
  end if
  
  vr_sign = vort500_max - vort200_max
  
  if (vr_sign .le. 0) then
    print *, 'X.X.X.X - VORT INCREASES BW 500 and 200: ',vr_sign
    print *, 'X.X.X.X - 500 ',vort500_max,'     200 ',vort200_max
    CYCLE
  end if
  
  PRINT *, 'O.O.O.O - SATISFIED VORT SLOPE CRITERIA'  

!-------------------------------------------------------------------
! --- STEP 4: LOCATE LOCAL SEA LEVEL PRESSURE MIN 
!-------------------------------------------------------------------
  
! Set large number for initial psmn
  psmn = 1E7
  psmin_i = 1
  psmin_j = 1
 
! Find starting indices for loop (top left corner of square)
  starti = i-psl_crit_dist
  startj = j-psl_crit_dist
  
! Loop over gridboxes within critdist to see where psl_min is
  do ii = i-psl_crit_dist,i+psl_crit_dist
  
  !write (*,'(<33>(I1,X))') ( pslVortColMtx(ii-starti+1,jj-startj+1), jj=j-psl_crit_dist,j+psl_crit_dist )
  
  do jj = j-psl_crit_dist,j+psl_crit_dist
! Makes sure it's within circular domain defined by pslVortColMtx
  !print *,ii-starti+1,' ',jj-startj+1
  if (pslVortColMtx(ii-starti+1,jj-startj+1) .gt. 0) then
      if (psl(ii,jj).lt.psmn) then
        psmn = psl(ii,jj)
        psmin_i = ii
        psmin_j = jj
      else
      end if
  end if
  end do
  end do
  
  lon_psl = rlon(psmin_i)
  lat_psl = rlat(psmin_j)
  print *, '  FOUND PSL MIN at ',lon_psl,' ',lat_psl
  
  psl_min = psmn

!
! SAVE INFO
!

    number = number + 1
    
    print *, 'ooooooooooooo -- CANDIDATE CYCLONE'
    print *, number, lon_psl, lat_psl
    print *, ' '

  if( number > nsmax ) then
    PRINT *, '***************************************'
    PRINT *, '  GOT TOO MANY STORMS - INCREASE nsmax '
    PRINT *, '***************************************'
    STOP
  endif

        idex(number) = i - nx
        jdex(number) = j - nx

   svort_max(number) =  vort_max
   svort_lon(number) =  lon_vort
   svort_lat(number) =  lat_vort

   swind_max(number) =  wind_max

    spsl_min(number) =  psl_min
    spsl_lon(number) =  lon_psl
    spsl_lat(number) =  lat_psl

!  WRITE(*,*) i-nx, j-nx
! zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
  end do
  end do

!===================================================================
! ---  OUTPUT
!===================================================================

  WRITE(iucy,*) iday,  imon, iyear, number, ihour
  WRITE(   *,*) iyear, imon, iday,  ihour,  number

  if( number == 0 ) RETURN

! ---  sthick_is(1:number) = .true.

  do i = 1,number
    WRITE(iucy,*) idex(i),       jdex(i),              & 
              spsl_lon(i),   spsl_lat(i),              &
              swind_max(i), svort_max(i), spsl_min(i)
!                stwc_is(i), sthick_is(i)
!                stwc_is(i), sthick_is(i), spsl_fail(i), spsl_drop(i)
  end do
  
!  deallocate(twcPSLColMtx)
  deallocate(pslVortColMtx)
!  deallocate(ringmtx)
  
!===================================================================
  end SUBROUTINE TSTORMS_HODGES

!######################################################################

  SUBROUTINE TSTORMS_WINTER ( Gwind, Gvort, Gpsl, &
                       Grlon, Grlat, iyear, imon, iday, ihour, &
                       iucy )
                       
! This subroutine was added to deal with high resolution 0.25 deg
! or finer data that the spline routine takes too long to deal with
! Colin Zarzycki

!===================================================================
! --- LOCATE TROPICAL STORMS 
!===================================================================

  use TSGPAD_MOD, only : GPAD2,  GPAD1FLIP
  use SPLINE_MOD, only : SPLIE2, SPLIE3, FRPRM, SHAPE
  implicit none

!-------------------------------------------------------------------
! --- INPUT ARGUMENTS 
!     Gwind  - wind speed at 850 mb
!     Gvort  - vorticity  at 850 mb
!     Gpsl   - sea level pressure
!     Grlon  - longitudes
!     Grlat  - latitudes
!     iyear  - year
!     imon   - month
!     iday   - day of month
!     ihour  - hour
!     iucy   - unit for output
!-------------------------------------------------------------------
! --- OUTPUT - file "cyclones" 
!-------------------------------------------------------------------
! --- record # 1 
!     num0   - day
!     imon0  - month
!     iyear  - year
!     number - number of cyclones found
! --- records # 2...number+1 
!     idex, jdex - (i,j) index of cyclone 
!     svort_max  - max vorticity                  
!     swind_max  - max wind              
!      spsl_min  - min sea level pressure                 
!     svort_lon,  svort_lat - longitude & latitude of max vorticity 
!      spsl_lon,   spsl_lat - longitude & latitude of min slp 
!      stwc_lon,   stwc_lat - longitude & latitude of warm core 
!    sthick_lon, sthick_lat - longitude & latitude of max thickness 
!-------------------------------------------------------------------

  real,    intent(in),    dimension(:,:) :: Gwind, Gvort,  Gpsl
  real,    intent(in),    dimension(:)   :: Grlon, Grlat
  integer, intent(in)                    :: ihour, iday,  imon, iyear, iucy

!-------------------------------------------------------------------
! --- LOCAL
!-------------------------------------------------------------------

  integer, parameter :: nx   = gnx0
  !integer, parameter :: nx   = 16
  integer, parameter :: nx2  = 2*nx
  integer, parameter :: nxp1 = nx + 1

  real,    parameter :: ftol  = 0.01   
  integer, parameter :: nsmax = 10000

  real, dimension(SIZE(Grlon)+nx2) :: rlon
  real, dimension(SIZE(Grlat)+nx2) :: rlat

  real, dimension(SIZE(Gwind,1)+nx2,SIZE(Gwind,2)+nx2) ::  &
        vort, wind, psl,    &
                    psl_dx, &
                    psl_dy

  real    :: vort_max, wind_max, psl_min
  real    :: lon_vort,           lon_psl
  real    :: lat_vort,           lat_psl

  integer, dimension(nsmax) :: idex, jdex

  real,    dimension(nsmax) ::                                    &
           svort_max, swind_max,  spsl_min, &
           svort_lon,             spsl_lon, &
           svort_lat,             spsl_lat, &
           spsl_fail

  real,    dimension(2)  :: p

  real :: xx, yy, rr, fret

  integer :: ierr_pos, ierr_mag

  integer :: i, im, ip, ix, ixp3, ixp6
  integer :: j, jm, jp, jx, jxp3, jxp6, jxp6h 
  integer :: number, iter
  
!====== Colin stuff

  integer :: starti,startj
  
  real :: distbwpts,avg

  integer :: ii,jj
  
  integer :: psmin_i,psmin_j,psl_crit_dist,num_fail_psl
  integer :: psl_off_dist, psl_od_diag
  real    :: psmn

!===================================================================

  ix    = SIZE( Gwind, 1 )
  jx    = SIZE( Gwind, 2 )
  ixp3  = ix + nx
  jxp3  = jx + nx
  ixp6  = ix + nx2
  jxp6  = jx + nx2
  jxp6h = jxp6 / 2

  number = 0

!-------------------------------------------------------------------
! --- SETUP
!-------------------------------------------------------------------

! Here we are "padding" the arrays with ghost cells so we can do the
! vorticity and wind checks at all locations

  CALL GPAD2( Gwind,  wind  )    !  Wind speed at 850 mb 
  CALL GPAD2( Gvort,  vort  )    !  Vorticity  at 850 mb 
  CALL GPAD2( Gpsl,   psl   )    !  Sea level pressure.
  CALL GPAD1FLIP( Grlon,  rlon, 0  )    !  Longitudes
  CALL GPAD1FLIP( Grlat,  rlat, 1  )    !  Latitudes  
 
! --- change sign of vorticity in southern hemisphere
  do j=1,jxp6
  if( rlat(j) < 0.0 ) then
    vort(:,j) = -1.0 * vort(:,j) 
  end if
  end do

!===================================================================
! -- Define areas
!===================================================================
  
  print *, '---------------------'
  print *, 'resolution ', resolution
  print *, 'pslVortMaxDist ',pslVortMaxDist
  print *, 'crit_dist ', crit_dist
  print *, 'crit_vort ', crit_vort
  print *, 'crit_twc ', crit_twc
  print *, 'dist_twc ', dist_twc
  print *, 'crit_thick ', crit_thick
  print *, 'dist_thick ', dist_thick
  print *, 'crit_psl ', crit_psl
  print *, 'dist_psl ', dist_psl
  print *, '---------------------'
  
  psl_crit_dist = CEILING(pslVortMaxDist/resolution)
  print *, 'psl_crit_dist ', psl_crit_dist

  psl_off_dist = CEILING(dist_psl/resolution)
  psl_od_diag = CEILING(psl_off_dist/SQRT(2.))
  print *, 'psl_off_dist ', psl_off_dist
  print *, 'psl_od_diag ', psl_od_diag
  
!===================================================================
! -- LOOP OVER GRID & LOOK FOR STORMS
!===================================================================
  
  do j = nxp1,jxp3
         if( ( rlat(j) > lat_bound_n ) .or. &
             ( rlat(j) < lat_bound_s ) ) CYCLE
  do i = nxp1,ixp3
! zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz

    im = i - nx
    ip = i + nx
    jm = j - nx
    jp = j + nx

!-------------------------------------------------------------------
! --- STEP 1: CHECK FOR VORTICITY MAX
!-------------------------------------------------------------------

   vort_max = MAXVAL( vort(im:ip,jm:jp) )
   wind_max = MAXVAL( wind(im:ip,jm:jp) )

! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  if( ( vort(i,j) /= vort_max  )  .or. &
      ( vort(i,j) <  crit_vort ) ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

  lon_vort = rlon(i)
  lat_vort = rlat(j)
  
  PRINT *, '        '
  PRINT *, '----------------------'
  PRINT *, 'O.O.O.O - FOUND VORT MAX at ',lon_vort,' ',lat_vort

!-------------------------------------------------------------------
! --- STEP 2: LOCATE LOCAL SEA LEVEL PRESSURE MIN 
!-------------------------------------------------------------------

  ierr_pos  = 0
  
! Set large number for initial psmn
  psmn = 1E7
  psmin_i = 1
  psmin_j = 1
 
! Find starting indices for loop (top left corner of square)
  starti = i-psl_crit_dist
  startj = j-psl_crit_dist

  call find_minmax_gen (.true.,psl(i-psl_crit_dist:i+psl_crit_dist,j-psl_crit_dist:j+psl_crit_dist),ii,jj,avg)
  psmin_i = ii+starti-1
  psmin_j = jj+startj-1
  psmn = psl(psmin_i,psmin_j)

  lon_psl = rlon(psmin_i)
  lat_psl = rlat(psmin_j)
  print *, '  FOUND PSL MIN at ',lon_psl,' ',lat_psl

  ! Calculate distance between vort center (i,j) and psl (psmin_i,psmin_j)
  ! in grid boxes
  
  !distbwpts = SQRT( (psmin_i-i)**2. + (psmin_j-j)**2. )
  !if (distbwpts .gt. psl_crit_dist) then
  !  print *, 'X.X.X.X - PSL dist = ',distbwpts,' compared to ',psl_crit_dist
  !  ierr_pos = 1
  !else
  !  print *, 'O.O.O.O - PSL dist = ',distbwpts,' compared to ',psl_crit_dist
  !  ierr_pos = 0
  !end if

  distbwpts = DGCDIST(lat_psl,lon_psl,lat_vort,lon_vort,2)
  if (distbwpts .gt. pslVortMaxDist) then
    print *, 'X.X.X.X - PSL dist = ',distbwpts,' compared to ',pslVortMaxDist
    ierr_pos = 1
  else
    print *, 'O.O.O.O - PSL dist = ',distbwpts,' compared to ',pslVortMaxDist
    ierr_pos = 0
  end if
  
  psl_min = psmn
  
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
   if ( ierr_pos == 1 ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

  num_fail_psl=0
  ierr_mag=0
  call check_contour(.true.,psl,psmin_i,psmin_j,num_fail_psl,1,ierr_mag,psl_off_dist,psl_od_diag,crit_psl)
  
  if ( ierr_mag == 1 ) then
    print *, 'X.X.X.X - PSL NUMFAIL = ',num_fail_psl,' compared to 0'
  else
    print *, 'O.O.O.O - PSL NUMFAIL = ',num_fail_psl,' compared to 0'
  end if

! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
   if ( ierr_mag == 1 ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

!-------------------------------------------------------------------
! --- WE HAVE A STORM.
!-------------------------------------------------------------------

    number = number + 1
    
    print *, 'ooooooooooooo -- CANDIDATE CYCLONE'
    print *, number, lon_psl, lat_psl!, '     ', lon_twc, lat_twc
    print *, ' '

  if( number > nsmax ) then
    PRINT *, '***************************************'
    PRINT *, '  GOT TOO MANY STORMS - INCREASE nsmax '
    PRINT *, '***************************************'
    STOP
  endif

        idex(number) = i - nx
        jdex(number) = j - nx

   svort_max(number) =  vort_max
   svort_lon(number) =  lon_vort
   svort_lat(number) =  lat_vort

   swind_max(number) =  wind_max

    spsl_min(number) =  psl_min
    spsl_lon(number) =  lon_psl
    spsl_lat(number) =  lat_psl

  end do
  end do

!===================================================================
! ---  OUTPUT
!===================================================================

  WRITE(iucy,*) iday,  imon, iyear, number, ihour
  WRITE(   *,*) iyear, imon, iday,  ihour,  number

  if( number == 0 ) RETURN

  do i = 1,number
    WRITE(iucy,100) idex(i), jdex(i), spsl_lon(i), spsl_lat(i), &
              swind_max(i), svort_max(i), spsl_min(i)

100  FORMAT(I6,1x,I6,1x,2(f7.2,2x),F7.3,2x,ES10.3,2x,F10.2)

  end do
  
!===================================================================
  end SUBROUTINE TSTORMS_WINTER

  SUBROUTINE TSTORMS_WINTER_POLE ( Gwind, Gvort, Gpsl, &
                       Grlon, Grlat, iyear, imon, iday, ihour, &
                       iucy )
                       
! This subroutine was added to deal with high resolution 0.25 deg
! or finer data that the spline routine takes too long to deal with
! Colin Zarzycki

!===================================================================
! --- LOCATE TROPICAL STORMS 
!===================================================================

  use TSGPAD_MOD, only : GPAD2_VAR,  GPAD1
  use SPLINE_MOD, only : SPLIE2, SPLIE3, FRPRM, SHAPE
  implicit none

!-------------------------------------------------------------------
! --- INPUT ARGUMENTS 
!     Gwind  - wind speed at 850 mb
!     Gvort  - vorticity  at 850 mb
!     Gpsl   - sea level pressure
!     Grlon  - longitudes
!     Grlat  - latitudes
!     iyear  - year
!     imon   - month
!     iday   - day of month
!     ihour  - hour
!     iucy   - unit for output
!-------------------------------------------------------------------
! --- OUTPUT - file "cyclones" 
!-------------------------------------------------------------------
! --- record # 1 
!     num0   - day
!     imon0  - month
!     iyear  - year
!     number - number of cyclones found
! --- records # 2...number+1 
!     idex, jdex - (i,j) index of cyclone 
!     svort_max  - max vorticity                  
!     swind_max  - max wind              
!      spsl_min  - min sea level pressure                 
!     svort_lon,  svort_lat - longitude & latitude of max vorticity 
!      spsl_lon,   spsl_lat - longitude & latitude of min slp 
!      stwc_lon,   stwc_lat - longitude & latitude of warm core 
!    sthick_lon, sthick_lat - longitude & latitude of max thickness 
!-------------------------------------------------------------------

  real,    intent(in),    dimension(:,:) :: Gwind, Gvort,  Gpsl
  real,    intent(in),    dimension(:)   :: Grlon, Grlat
  integer, intent(in)                    :: ihour, iday,  imon, iyear, iucy

!-------------------------------------------------------------------
! --- LOCAL
!-------------------------------------------------------------------

  integer, parameter :: nx   = gnx0
  !integer, parameter :: nx   = 16
  integer, parameter :: nx2  = 2*nx
  integer, parameter :: nxp1 = nx + 1

  real,    parameter :: ftol  = 0.01   
  integer, parameter :: nsmax = 10000

  integer :: nlon, nlat, nlond2, varnx, iloopst

  real, dimension(SIZE(Grlon)*3) :: rlon
  real, dimension(SIZE(Grlat)+nx2) :: rlat

  real, dimension(SIZE(Gwind,1)*3,SIZE(Gwind,2)+nx2) ::  &
        vort, wind, psl,    &
                    psl_dx, &
                    psl_dy

  real    :: vort_max, wind_max, psl_min
  real    :: lon_vort,           lon_psl
  real    :: lat_vort,           lat_psl

  integer, dimension(nsmax) :: idex, jdex

  real,    dimension(nsmax) ::                                    &
           svort_max, swind_max,  spsl_min, &
           svort_lon,             spsl_lon, &
           svort_lat,             spsl_lat, &
           spsl_fail
           
  logical, dimension(nsmax) ::                                    &
                                            stwc_is,  sthick_is

  real,    dimension(2)  :: p

  real :: xx, yy, rr, fret

  integer :: ierr_pos, ierr_mag

  integer :: i, im, ip, ix, ixp3, ixp6
  integer :: j, jm, jp, jx, jxp3, jxp6, jxp6h 
  integer :: number, iter
  
!====== Colin stuff

  integer :: starti,startj
  
  real :: distbwpts, avg, psl_grad

  integer :: ii,jj
  
  integer :: psmin_i,psmin_j,psl_crit_dist,num_fail_psl
  integer :: psl_off_dist, psl_od_diag
  real    :: psmn

!===================================================================

  ix    = SIZE( Gwind, 1 )
  jx    = SIZE( Gwind, 2 )
  ixp3  = ix + SIZE( Gwind, 1 )
  jxp3  = jx + nx
  ixp6  = ix + SIZE( Gwind, 1 )*2
  jxp6  = jx + nx2
  jxp6h = jxp6 / 2

  nlat = SIZE( Gwind, 2 )
  nlon = SIZE( Gwind, 1 )
  nlond2 = SIZE( Gwind, 1 ) / 2

  number = 0

!-------------------------------------------------------------------
! --- SETUP
!-------------------------------------------------------------------

! Here we are "padding" the arrays with ghost cells so we can do the
! vorticity and wind checks at all locations

  CALL GPAD2_VAR( Gwind,  wind  )    !  Wind speed at 850 mb 
  CALL GPAD2_VAR( Gvort,  vort  )    !  Vorticity  at 850 mb 
  CALL GPAD2_VAR( Gpsl,   psl   )    !  Sea level pressure.
  CALL GPAD1( Grlon,  rlon  )    !  Longitudes
  CALL GPAD1( Grlat,  rlat  )    !  Latitudes  

! --- change sign of vorticity in southern hemisphere
  do j=1,jxp6
  if( rlat(j) < 0.0 ) then
    vort(:,j) = -1.0 * vort(:,j) 
  end if
  end do

!===================================================================
! -- Define areas
!===================================================================
  
  print *, '---------------------'
  print *, 'resolution ', resolution
  print *, 'pslVortMaxDist ',pslVortMaxDist
  print *, 'crit_dist ', crit_dist
  print *, 'crit_vort ', crit_vort
  print *, 'crit_twc ', crit_twc
  print *, 'dist_twc ', dist_twc
  print *, 'crit_thick ', crit_thick
  print *, 'dist_thick ', dist_thick
  print *, 'crit_psl ', crit_psl
  print *, 'dist_psl ', dist_psl
  print *, '---------------------'

      psl_crit_dist = CEILING(pslVortMaxDist/resolution)
      print *, 'psl_crit_dist ', psl_crit_dist

      psl_off_dist = CEILING(dist_psl/resolution)
      psl_od_diag = CEILING(psl_off_dist/SQRT(2.))
      print *, 'psl_off_dist ', psl_off_dist
      print *, 'psl_od_diag ', psl_od_diag
  

  
!===================================================================
! -- LOOP OVER GRID & LOOK FOR STORMS
!===================================================================

  print *, 'nx ', nx
  print *, 'nx2 ', nx2
  print *, 'nxp1 ', nxp1
  print *, 'ix ', ix
  print *, 'jx ', jx
  print *, 'ixp3 ', ixp3
  print *, 'jxp3 ', jxp3
  print *, 'ixp6 ', ixp6
  print *, 'jxp6 ', jxp6
  print *, 'jxp6h ', jxp6h

  iloopst = ix + 1

  do j = nxp1,jxp3
         if( ( rlat(j) > lat_bound_n ) .or. &
             ( rlat(j) < lat_bound_s ) ) CYCLE

     varnx = SIND(ABS(rlat(j)))*(nlond2-nx) + nx

  do i = iloopst,ixp3
! zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz

    im = i - varnx
    ip = i + varnx
    jm = j - nx
    jp = j + nx

!-------------------------------------------------------------------
! --- STEP 1: CHECK FOR VORTICITY MAX
!-------------------------------------------------------------------

   vort_max = MAXVAL( vort(im:ip,jm:jp) )
   wind_max = MAXVAL( wind(im:ip,jm:jp) )

! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  if( ( vort(i,j) /= vort_max  )  .or. &
      ( vort(i,j) <  crit_vort ) ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

  lon_vort = rlon(i)
  lat_vort = rlat(j)
  
  PRINT *, '        '
  PRINT *, '----------------------'
  PRINT *, ' ',vort_max
  PRINT *, 'O.O.O.O - FOUND VORT MAX at ',lon_vort,' ',lat_vort

!-------------------------------------------------------------------
! --- STEP 2: LOCATE LOCAL SEA LEVEL PRESSURE MIN 
!-------------------------------------------------------------------

  ierr_pos  = 0
  
! Set large number for initial psmn
  psmn = 1E7
  psmin_i = 1
  psmin_j = 1

  call find_minmax_gen (.true., psl(im:ip,jm:jp), ii, jj, avg)
  print *,ip-im
  print *,jp-jm
  !print *,avg/100.
  psmin_i = ii+im-1
  psmin_j = jj+jm-1  
  psmn = psl(psmin_i,psmin_j)
  lon_psl = rlon(psmin_i)
  lat_psl = rlat(psmin_j)

  xx      = lon_psl - lon_vort
  yy      = lat_psl - lat_vort
  distbwpts      = SQRT(xx*xx + yy*yy)

  if (distbwpts .gt. pslVortMaxDist) then
    print *, 'X.X.X.X - PSL dist = ',distbwpts,' compared to ',pslVortMaxDist
    ierr_pos = 1
  else
    print *, 'O.O.O.O - PSL dist = ',distbwpts,' compared to ',pslVortMaxDist
    ierr_pos = 0
  end if

! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
   !if ( ierr_pos == 1 ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

  psl_min = psmn

  psl_grad = avg - psl_min

  if (psl_grad .lt. crit_psl) then
    print *, 'X.X.X.X - PSL grad = ',psl_grad,' compared to ',crit_psl
    ierr_pos = 1
  else
    print *, 'O.O.O.O - PSL grad = ',psl_grad,' compared to ',crit_psl
    ierr_pos = 0
  end if

! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
   if ( ierr_mag == 1 ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

!-------------------------------------------------------------------
! --- WE HAVE A STORM.
!-------------------------------------------------------------------

    number = number + 1
    
    print *, 'ooooooooooooo -- CANDIDATE CYCLONE'
    print *, number, lon_psl, lat_psl!, '     ', lon_twc, lat_twc
    print *, ' '

  if( number > nsmax ) then
    PRINT *, '***************************************'
    PRINT *, '  GOT TOO MANY STORMS - INCREASE nsmax '
    PRINT *, '***************************************'
    STOP
  endif

        idex(number) = i - nx
        jdex(number) = j - nx

   svort_max(number) =  vort_max
   svort_lon(number) =  lon_vort
   svort_lat(number) =  lat_vort

   swind_max(number) =  wind_max

    spsl_min(number) =  psl_min
    spsl_lon(number) =  lon_psl
    spsl_lat(number) =  lat_psl

  end do
  end do

!===================================================================
! ---  OUTPUT
!===================================================================

  WRITE(iucy,*) iday,  imon, iyear, number, ihour
  WRITE(   *,*) iyear, imon, iday,  ihour,  number

  if( number == 0 ) RETURN

  do i = 1,number
    WRITE(iucy,100) idex(i), jdex(i), spsl_lon(i), spsl_lat(i), &
              swind_max(i), svort_max(i), spsl_min(i)

100  FORMAT(I6,1x,I6,1x,2(f7.2,2x),F7.3,2x,ES10.3,2x,F10.2)

  end do
  
!===================================================================
  end SUBROUTINE TSTORMS_WINTER_POLE

!######################################################################
!######################################################################

  SUBROUTINE find_minmax (do_min, var, imm, jmm, avg)
    implicit none

    logical, intent(in)                    :: do_min
    real,    intent(in),    dimension(2*gnx0+1,2*gnx0+1) :: var
    integer, intent(out)                   :: imm, jmm
    real,    intent(out)                   :: avg

    real    ::  varminmax
    integer ::  i,j

    if (do_min) then
       varminmax = MINVAL( var(:,:) )
    else
       varminmax = MAXVAL( var(:,:) )
    end if
    avg = 0.
    do i = 1, 2*gnx0+1
       do j = 1, 2*gnx0+1
          avg = avg + var(i,j)
          if (var(i,j) == varminmax) then
             imm = i
             jmm = j
          end if
       end do
    end do
    avg = avg / ((2.*gnx0+1)*(2.*gnx0+1))
  end SUBROUTINE find_minmax

!######################################################################
  SUBROUTINE find_minmax_gen (do_min, var, imm, jmm, avg)
    implicit none

    logical, intent(in)                    :: do_min
    real,    intent(in),    dimension(:,:) :: var
    integer, intent(out)                   :: imm, jmm
    real, intent(out) :: avg

    real    ::  varminmax
    integer ::  i,j
    integer ::  nx,ny,ntot

    nx = SIZE(var,1)
    ny = SIZE(var,2)
    ntot = nx*ny

    avg = 0.
    if (do_min) then
       varminmax = MINVAL( var(:,:) )
    else
       varminmax = MAXVAL( var(:,:) )
    end if

    do i = 1, SIZE(var,1)
       do j = 1, SIZE(var,2)
          avg = avg + var(i,j)
          if (var(i,j) == varminmax) then
             imm = i
             jmm = j
          end if
       end do
    end do
    avg = avg / float(ntot)
    print *,avg
  end SUBROUTINE find_minmax_gen

!######################################################################

REAL FUNCTION DGCDIST(RLAT1,RLON1,RLAT2,RLON2,IU)
  IMPLICIT NONE


! calculate the great circle distance between two points
! usage: dist = gcdist (rlat1,rlon1,rlat2,rlon2,iu)
! nomenclature :
! .   rlat1,rlon1 - latitude and longtitude of the first point
! .   rlat2,rlon2 - latitude and longtitude of the second point
! .   iu          - code for the type units gcdist is to return
! .               = 1 : gcdist returned in radians
! .               = 2 : gcdist returned in degrees
! .               = 3 : gcdist returned in meters
! .               = 4 : gcdist returned in kilometers
! .               = 5 : gcdist returned in *not used*

  integer, intent(in) :: IU
  real, intent(in) :: RLAT1, RLON1, RLAT2, RLON2

  real :: UNITS(5), RAD, DLONR, RLAT1R, RLAT2R
! local stuff
  DATA UNITS/1.0D0,57.29577995691645D0,6371220.D0,6371.2200D0,0.D0/
! change as required
  DATA RAD/0.01745329238474369D0/

! special test if RLAT1=RLAT2 and RLON1=RLON2
      IF(RLAT1.EQ.RLAT2.AND.RLON1.EQ.RLON2) THEN
         DGCDIST = 0.D0
         RETURN
      END IF
      RLAT1R = RLAT1*RAD
      RLAT2R = RLAT2*RAD
      DLONR = DMIN1(ABS(RLON1-RLON2),ABS(360.D0-RLON1+RLON2) &
     +        ABS(360.D0-RLON2+RLON1))*RAD

      DGCDIST = ACOS(SIN(RLAT1R)*SIN(RLAT2R) &
     +          COS(RLAT1R)*COS(RLAT2R)*COS(DLONR))*UNITS(IU)

end function DGCDIST


!###################################################################################
  SUBROUTINE check_contour(do_min,var_in,mxmn_i,mxmn_j,num_fail,max_fail,ierr,off_dist,od_diag,crit_df)

  ! This subroutine takes a 2-D field and checks if there is a closed contour
  ! of sufficient critical falloff (or fallup) in up to 8 quadrants around a max/min
  ! Colin Zarzycki 10/4/2014

  ! do_min (logical) -- .false. if using max, .true. if using min
  ! var_in -- full 2-D field
  ! mxmin_i, mxmin_j -- absolute i,j index of FULL field (from 0,0)
  ! num_fail -- number of quadrants that fail to have closed contour
  ! max_fail -- max allowable quadrants to fail (0 means fully closed contour must be present)
  ! ierr -- returns 0 if num_fail/max_fail satisfied, 1 if failed
  ! off_dist -- max lateral allowable distance in absolute cells
  ! od_diag  -- max diagonal allowable distance in absolute cells
  ! crit_dif -- difference of contour from max (i.e., 4 mb in PSL)

    implicit none

    !real,    intent(in),    dimension(2*gnx0+1,2*gnx0+1) :: var
    !real, intent(in), dimension(744,385) :: var_in
    !ix+nx2, jx+nx2
    real,    intent(in), dimension(:,:)    :: var_in
    integer, intent(in)                    :: mxmn_i, mxmn_j
    integer, intent(out)                   :: num_fail
    integer, intent(out)                   :: ierr
    integer, intent(in)                    :: off_dist, od_diag
    real,    intent(in)                    :: crit_df
    logical, intent(in)                    :: do_min
    integer, intent(in)                    :: max_fail

    integer :: mm, mmm
    real :: vr_n, vr_s, vr_e, vr_w
    real :: vr_nw, vr_sw, vr_ne, vr_se
    real :: vr_n_temp, vr_s_temp, vr_e_temp, vr_w_temp
    real :: vr_nw_temp, vr_sw_temp, vr_ne_temp, vr_se_temp

    real, dimension(SIZE(var_in,1),SIZE(var_in,2)) :: var

    if (do_min) then
       var = -var_in
    else
       var = var_in
    end if

! This is a first guess at falloff by looking at max acceptable radius

  vr_n = var(mxmn_i,mxmn_j) - var(mxmn_i - off_dist, mxmn_j)
  vr_s = var(mxmn_i,mxmn_j) - var(mxmn_i + off_dist, mxmn_j)
  vr_e = var(mxmn_i,mxmn_j) - var(mxmn_i, mxmn_j - off_dist)
  vr_w = var(mxmn_i,mxmn_j) - var(mxmn_i, mxmn_j + off_dist)
  !! Cross terms
  vr_nw = var(mxmn_i,mxmn_j)-var(mxmn_i-od_diag,mxmn_j-od_diag)
  vr_sw = var(mxmn_i,mxmn_j)-var(mxmn_i+od_diag,mxmn_j-od_diag)
  vr_ne = var(mxmn_i,mxmn_j)-var(mxmn_i-od_diag,mxmn_j+od_diag)
  vr_se = var(mxmn_i,mxmn_j)-var(mxmn_i+od_diag,mxmn_j+od_diag)

! Loop over radii inside of max acceptable radius to see if there are any other values
! which are closed contours

  do mm = 1,off_dist
  vr_n_temp = var(mxmn_i,mxmn_j) - var(mxmn_i - mm, mxmn_j)
  vr_s_temp = var(mxmn_i,mxmn_j) - var(mxmn_i + mm, mxmn_j)
  vr_e_temp = var(mxmn_i,mxmn_j) - var(mxmn_i, mxmn_j - mm)
  vr_w_temp = var(mxmn_i,mxmn_j) - var(mxmn_i, mxmn_j + mm)
  if (vr_n_temp.gt.vr_n) then
  	vr_n = vr_n_temp
  end if
  if (vr_s_temp.gt.vr_s) then
  	vr_s = vr_s_temp
  end if
  if (vr_e_temp.gt.vr_e) then
  	vr_e = vr_e_temp
  end if
  if (vr_w_temp.gt.vr_w) then
  	vr_w = vr_w_temp
  end if
  end do

  do mmm = 1,od_diag
  vr_nw_temp = var(mxmn_i,mxmn_j) - var(mxmn_i - mmm, mxmn_j - mmm)
  vr_sw_temp = var(mxmn_i,mxmn_j) - var(mxmn_i + mmm, mxmn_j - mmm)
  vr_ne_temp = var(mxmn_i,mxmn_j) - var(mxmn_i - mmm, mxmn_j + mmm)
  vr_se_temp = var(mxmn_i,mxmn_j) - var(mxmn_i + mmm, mxmn_j + mmm)
  if (vr_nw_temp.gt.vr_nw) then
  	vr_nw = vr_nw_temp
  end if
  if (vr_sw_temp.gt.vr_sw) then
  	vr_sw = vr_sw_temp
  end if
  if (vr_ne_temp.gt.vr_ne) then
  	vr_ne = vr_ne_temp
  end if
  if (vr_se_temp.gt.vr_se) then
  	vr_se = vr_se_temp
  end if
  end do

  num_fail = 0

  ! Check if quadrants do not meet criteria

  if (vr_n.lt.crit_df) then
  	num_fail = num_fail + 1
  end if
  if (vr_s.lt.crit_df) then
  	num_fail = num_fail + 1
  end if
  if (vr_e.lt.crit_df) then
  	num_fail = num_fail + 1
  end if
  if (vr_w.lt.crit_df) then
  	num_fail = num_fail + 1
  end if
  if (vr_nw.lt.crit_df) then
  	num_fail = num_fail + 1
  end if
  if (vr_sw.lt.crit_df) then
  	num_fail = num_fail + 1
  end if
  if (vr_ne.lt.crit_df) then
  	num_fail = num_fail + 1
  end if
  if (vr_se.lt.crit_df) then
  	num_fail = num_fail + 1
  end if
  
  if (num_fail.gt.max_fail) then
  	ierr = 1
  else
    ierr = 0
  end if

  end SUBROUTINE check_contour

!#################################################################################

  SUBROUTINE circleMatrix (matrix, thresh_grid_dist)
    implicit none

    integer, intent(out), dimension(:,:)   :: matrix
    integer, intent(in)                    :: thresh_grid_dist

    integer :: iii,jjj,center_i,center_j
    real    :: dist
    real, parameter :: eps=10e-15

    center_i = thresh_grid_dist+1
    center_j = thresh_grid_dist+1
  
  do iii = 1,2*thresh_grid_dist+1
    do jjj = 1,2*thresh_grid_dist+1
      dist = ( (iii-center_i)**2 + (jjj-center_j)**2 )
      dist = dist ** (0.5)
      if (dist .LE. (thresh_grid_dist + eps)) then
        matrix(iii,jjj) = 1
      else
        matrix(iii,jjj) = 0
      end if
    end do
  end do

!  do iii=1,2*thresh_grid_dist+1
!    write (*,'(<2*thresh_grid_dist+1>(I1,X))') ( pslVortColMtx(iii,jjj), jjj=1,2*thresh_grid_dist+1 )
!  enddo

  end SUBROUTINE circleMatrix

!#################################################################################

!--  SUBROUTINE ringMatrix (ringmtx, matrix, thresh_grid_dist)
!--    implicit none
!--
!--    integer, intent(out), dimension(:,:)   :: ringmtx
!--    integer, intent(in), dimension(:,:)    :: matrix
!--    integer, intent(in)                    :: thresh_grid_dist
!--
!--    integer :: iii,jjj,count
!--    real    :: dist
!--    real, parameter :: eps=10e-15
!--
!--!  allocate ( ringmtx(2*thresh_grid_dist+1,2*thresh_grid_dist+1) )
!--!  ringmtx = pslVortColMtx
!--!
!--!  
!--!    do iii = 1,2*thresh_grid_dist+1
!--!        do jjj = 1,2*thresh_grid_dist+1
!--!            count = 1
!--!            if (matrix(iii,jjj) .EQ. 0) then
!--!                count = count+1
!--!            else
!--!                ringmtx(iii,jjj) = 2
!--!                ringmtx(iii,2*thresh_grid_dist+1-jjj+1) = 2
!--!                exit
!--!            end if
!--!        end do
!--!    end do
!--!    
!--!    do jjj = 1,2*thresh_grid_dist+1
!--!        do iii = 1,2*thresh_grid_dist+1
!--!            count = 1
!--!            if (matrix(iii,jjj) .EQ. 0) then
!--!                count = count+1
!--!            else if (matrix(iii,jjj) .EQ. 1) then
!--!                ringmtx(iii,jjj) = 2
!--!                ringmtx(2*thresh_grid_dist+1-iii+1,jjj) = 2
!--!                exit
!--!            else
!--!            end if
!--!        end do
!--!    end do
!--!
!--!  do iii=1,2*thresh_grid_dist+1
!--!    print *,( ringmtx(iii,jjj), jjj=1,2*thresh_grid_dist+1 )
!--!  enddo
!--
!--  end SUBROUTINE ringMatrix

!-------------------------------------------------------------------
! --- CREATE RING MATRIX (COLIN)
!-------------------------------------------------------------------

  
!#################################################################################

  end MODULE TSTORMS_MOD
